A Method for Calculating the Structure of (Singular) Spacetimes in the Large 

A formalism and its numerical implementation is presented which allows to calculate quantities 
determining the spacetime structure in the large directly. This is achieved by conformal techniques 
by which future null infinity {J^) and future timelike infinity (?"*") are mapped to grid points on 
the numerical grid. The determination of the causal structure of singularities, the localization of 
event horizons, the extraction of radiation, and the avoidance of unphysical reflections at the outer 
boundary of the grid, are demonstrated with calculations of spherically symmetric models with a 
scalar field as matter and radiation model. 
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I. INTRODUCTION 
A. General framework 
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From the singularity theorems by S. Hawking and R. Penrose it is known that the appearance of singularities in 
general relativity is an unavoidable feature for "strong" initial data jl| . Conversely, it has been proven recently that 
for small initial data the future of the initial value hypersurface looks in the large in a certain sense like the future of 
flat space data |^|-^. Unfortunately, it has turned out to be extremely difficult to obtain reasonable conjectures — 
^ , not to mention the proofs — about the properties of a spacetime given by an initial value problem, as illustrated, for 
example, by the history of the cosmic censorship hypothesis. 
\ Getting an overview for the phenomena appearing and making reasonable conjectures is an area numerical relativity 
^ 1 can contribute a lot to and actually already has; most impressive illustrated by the discovery of the echoing property 
Q^ by M. Choptuik |^. But most of the codes used so far are designed to analyze local behaviour. Statements about the 
CN , behaviour in the large, global issues, are obtained by assuming that the extension of the grid used approximates an 
infinite grid well enough. This way it is very difficult to get a reliable error estimate and to decide, what is sufficiently 
far, if not impossible at all. 

But what global issues are of interest and why are they interesting? 
\ Firstly, there are questions related to cosmic censorship, especially the location and penetration of an event horizon, 
^ ' which is the boundary of the region of spacetime from where no null geodesic reaches null infinity. In a recent review 
article about cosmic censorship C. Clarke writes about the location of event horizons: "In terms of numerical 
I simulations, this means that it is essential to perform the simulations in the compactified picture in which J'^ is 
represented explicitly" . 

Secondly, there is the whole topic of gravitational radiation on asymptotically fiat spacetimes. Due to the gauge 
^ ' freedom and the non-linearity of the theory the classification of radiation as ingoing and outgoing is a difficult issue. 
k>( ' In- and outgoing gravitational waves are defined with respect to null infinity only. The related difficulties in extracting 
' gravitational waves from the grid at finite distances and the problem with the avoidance of an unphysical refiection 
^ on the outer boundary of the grid are well 

known and have been a topic for research for a long time. In most methods developed so far the error made consists 
of the error from reading off at finite distances and the discretization error. 

In this paper I will present the numerical implementation of a formalism allowing to calculate a "compactified" 
spacetime including Sf^ and The solution of the problems concerning radiation is trivial by construction as 
will be seen. The only errors which may appear are caused by the discretization error of the numerical partial 
differential equation solver. These errors can be estimated and controlled by grid extrapolation techniques (e. g. 
Richardson extrapolation). Although the calculations have been done under the assumption of spherical symmetry 
with a scalar field as model for radiation the simplicity and exactness of radiation extraction hold for arbitrary 
symmetry assumptions. 

Furthermore, due to the spherical symmetry a special coordinate gauge could be found allowing to cover the complete 
domain of dependence of the initial hypersurface and to calculate the causal structure of a singularity. The location 
of the event horizon is straightforward. 



*This work is to a large extent part of my Ph. D. thesis which has been done at the Max-Planck-Institut fiir Astrophysik, 
Postfach 1523, D-85740 Garching 
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The formalism is based on conformal techniques developed by R. Penrose to describe asymptotically flat spacetimes. 
By a rescaling gab = dab the physical spacetime (M,^^^) is mapped to an unphysical spacetime {M, gab,^) with 
boundary J'. The boundary represents null infinity. M is a "compactifLed"[| version of M. Gravitational radiation is 
the value of certain components of the Weyl spinor on that boundary. In the originally suggested form the formalism 
is not suited to describe initial value problems. H. Friedrich has modified it to describe initial value problems — one 
has to solve a set of evolution equations for the unphysical spacetime {M, gab, The rescaling factor fl acts like an 
artificial matter field for the Einstein equations. Subsection |IIE| reviews the equations for the unphysical spacetime. 
J. Winicour and R. Gomez have developed an approach which uses Bondi's ideas for describing gravitational radiation: 
Outgoing null cones, with the area distance and the direction angle as coordinates on it, are compactified to represent 
future null infinity by grid points. Their method gives simpler equations but has certain disadvantages: The existence 
of smooth outgoing null cones is essential. Null caustics, which are caused e.g. by gravitational lensing effects, must 
not appear. Furthermore the evolution scheme cannot penetrate event horizons. The extraction of radiation is 
complicated slightly as the Weyl spinor is not a variable of the system — the determination of gravitational radiation 
requires to calculate numerically derivatives on J^. Their formalism solves a characteristic initial value, the formalism 
presented here solves a "normal" , spacelike initial value problems in unphysical spacetime. 

B. Description of asymptotically flat spacetimes with conformal techniques 

Shortly after H. Bondi et. al. Q proved that gravitational radiation is not a gauge effect R. Penrose suggested 
a coordinate independent way to characterize asymptotically flat spaces and gravitational radiation. A thorough 
discussion of the ideas and the interpretation can be found at various places in the literature, e.g. |^,^. The definitions 
of asymptotical flatness given in the literature differ slightly. The following will be used here: 

Definition 1 A spacetime {M,gab) is caHed asymptotically flat ij there is another "unphysical" spacetime {M,gab) 
with boundary J and a smooth embedding by which M can be identified with M — such that: 

1. There is a smooth function on M with 

^^Ia7>0 "'^d, gab i^'^gab- 

2. On J 

n = and Va^J ^ 0. 

3. Each null geodesic in {M,gab) acquires a past and a future endpoint on J . 

Because of item || null geodesically incomplete spacetimes like Schwarzschild are not asymptotically flat. The next 
definition includes those spacetimes which have only an asymptotically flat part: 

Definition 2 A spacetime is called weakly asymptotically flat if definition [J with the exception of item |^ is 
fulfilled. 

Definition |^ and |^ classify spacetimes, they do not require that the Einstein equation is fulfilled. For a physical 

problem one would like to give "asymptotically flat data" and evolve them according to the Einstein equation. 

Nevertheless the geometrically description was extremely helpful in analyzing asymptotically flat spacetimes and it 

has been successfully used as guideline to construct a formalism for analyzing initial value problems. This method 

has been developed and applied to various matter sources by H. Friedrich [p|jic|-p^ and myself 

The idea is to choose a spacelike initial value surface in the unphysical spacetime {M, gab) and to evolve it. 

For Minkowski space the unphysical spacetime {M,gab) can be smoothly extended with three points, future (i"*") and 

past timelike infinity, the end and the starting point of all timelike geodesies of (M, gab) respectively, and spacelike 

infinity (i"), the end point of all spacelike geodesies of (M, gab)- The point i° divides J into two disjunct parts, future 

[J^) and past {J~) null infinity. It is well known and has been discussed elsewhere that there are unsolved problems 



is not really compact — in Minkowski spacetime there are three points, future and past timclikc and spacelike infinity, 
missing. In general M cannot be smoothly extended to contain those points. 
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in smoothly extending a "normal" Cauchy hypersurface of M to if the spacetime has non-vanishing ADM mass. 
Certain curvature quantities blow up at ^ reflecting the non-invariance of the mass under rescalings. 
By choosing a spacelike (with respect to gab) hypersurface S not intersecting «° but {J ) we avoid the problems 
with . S is called a hyperboloidal hypersurface — the corresponding initial value problem is called a hyperboloidal 
initial value problem. 

The well-posedness of the hyperboloidal initial value problem for general relativistic scalar fields has been discussed 
in UJ. There, a precise definition of a hyperboloidal initial value problem can be found, too. 

Figure |l| shows a diagram of the unphysical spacetime with an example of a hyperboloidal surface in it. In figure |^ 
the corresponding physical spacetime is shown. In both figures only one space coordinate is drawn. All points, except 
those on the axis and represent spheres. The domain of dependence D{S) of S will not contain the complete 
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FIG. 1. Unphysical Minkowski spacetime 




spacetime. The interior of S corresponds to an everywhere spacelike hypersurface in the physical spacetime which 
approaches a null hypersurface N asymptotically. If is a light cone L then the domain of dependence of S is L. 
Therefore the hyperboloidal initial value problem is well suited to describe the future (past) of data on the spacelike 
hypersurface S, e. g. a stellar object and the gravitational radiation caused by its future time evolution. It is not well 
suited to investigate the structure near 

Why is it advantageous to solve the initial value problem in the unphysical spacetime? The future of the data in 
physical spacetime (M,^^^) is completely determined by the data on the interior of S. Since the rescaling factor 
f2 is known, the calculation of quantities of physical spacetime from the unphysical quantities is merely algebra. In 
unphysical spacetime Sf'^ is on the grid. J^^ is an ingoing null cone starting at the 2-surface in S where $7 vanishes. 
By extending S a little bit beyond the intersection with in unphysical spacetime we do not change the physical 
spacetime (M,5^{,). But with a finite number of gridpoints and thus a finite computation time the whole future of 
the initial hypersurface is covered. And since is an ingoing null line the values at M and ^7"*" — describing the 
physics — do not depend on the values at the outer boundary of the extended S. The numerical dependence on the 
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outer boundary must converge to 0. 

Future null infinity, JT""*", which can be found by searching for the = contour, is on the grid. By an appropriate 
choice of gauge it can even be arranged that the position of J on the grid is known by analytical considerations. In 
the worst case, one has to interpolate between neighbouring grid points to find the radiation emitted. In the case 
of gravitational radiation, it is given by certain components of the Weyl tensor, which are variables in the system of 
equations for the unphysical spacetime. 

II. THE MODEL AND THE SYSTEM OF EQUATIONS IN UNPHYSICAL SPACETIME 

A. The model 

The choice of the model is always a compromise between generality and manageability. For this work it was 
necessary to have some model for radiation to demonstrate the simplicity of the radiation extraction in the formalism. 
On the other side the spacetime symmetry should be as high as possible to reduce the demand for computational 
resources and to simplify coordinate choice related questions. 

The spacetime is assumed to be spherically symmetric, the radiational degree of freedom is modeled by a conformally 
invariant scalar field. The equations are 

60- 1^ = (la) 
6 

4'^ " /-ft- V • -r/v • 2 ■ ■ 4 

where the second is equivalent to Gah = n Tat with 



il-^n4>')Rab^ U(V,0)(V60)-iK0VaVb<^-iK.g„b(V'0)(V,0) ) , (lb) 



fab = (Va0)(Vb^) - i VaV60 + ^ 4>' Rab " \ Qab (^{^' 4>) c4>) + J 4>' Rj ■ (Ic) 

The notation used is explained in appendix |^. 

The form invariance of ( p^ ) under the rescalings gab — 9ab ^^"^ ^ ~ ^ reason for not choosing the 

massless Klein-Gordon scalar field with its simpler equations in physical spacetime. In |^ it has been shown that the 
initial value problems for these matter models are equivalent for practical purposes. 

B. The equations for unphysical spacetime 

In Q it has been argued that the Einstein equation, if simply translated to the unphysical spacetime, becomes 
"singular" on J7 and thus is not suitable to calculate the unphysical spacetime {M,gab)- For a derivation of regular 
equations for the unphysical spacetime the reader is referred to |^ , here only the final first order set of equations is 
given. 

In this first order system, the variables are the components of the frame Cfc" with respect to the coordinates x^, ei/^, 
the frame components j—jk of the Ricci rotation coefficients 7°jfc, certain combinations of the components of the 

tracefree part Rab of the Ricci tensor and the regularized Weyl tensor dabc'^ '■= Cabc'^, the conformal factor il, the 
frame components fli of its gradient Qa, the trace of its second derivative to, the conformal scalar field (j), the frame 

components (pi of its gradient (f>a and combinations of its second derivatives (pab ='■ <pab + \9ab<pc'- The first order 
system, written in abstract index notation is: 

AAe%c:=T%c = (2a) 

A/'-y abc" '■= Rdm abc'^ ~ Ralg abc^ = (2b) 

1 2 

A/'rq6c V[Qi?f,]c + -^{"^ [aR) 9b\c - ^ddab/ + TUabc ~ ^ ^ m[a\d\'^ gb]c 

= (2c) 
2 

^fd abc ■= ^ddab/ - mabc + ^ '7l[a|d|'' 5b]c = (2d) 

AAn a := Vaf^ - f^a = (2e) 
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AAdSI ah ■■= "^a^b +]^nRab-LOgab-]^^^ Tab = (2f) 

N^ab-=Va^+\Rab^'' + ^R^a + ^^VaR-\n'' Tab 

= (2g) 

J^4,a-=^aCf>-cl)a=0 (2h) 

A/'d0 ab ■= ^a4>b - 4>ab - -J^'c" 9ab = (2i) 



AAn0 <^a'^ - I - (2j) 

A/'dD0 abc := V[a0b]c + ^ (0V[ai? + i?(/)[a)g6]c " ^ Ralg ab/ (l^d = (2k) 
A^DD^a :=V,0fc''-i(0Vai? + i?0a)=O, (21) 

6 

where T^^c is the torsion whose components are given by 
the tensor i?aig abc'' is an abbreviation for 



-kj 



Ra\g abed = ^ dabcd + gc[aRb]d " gd[aRb]c + T 5c[a 9b]d R, (3) 



6 

^^diff abc"^ has the components 

-Rdiff ijfe- = ej(^-ik) — eiij-jj^ — ^-iml—j_k + l-j_ml—ik 

I m / I m / / /T-im 

^ ; r;' ! mk ~r / / m/c / mk jij 

the "energy momentum tensor" in unphysical spacetime is 

Tab = <f>a 4lb - (fl <l>ab + ^<l>'^ Rab - ^ ffab (^c (p" + -^(f^ R 

and 

1 

mabc 



1 - i 17202 



3 1 1 1 » 3 - 

[2 <^[a06]c - gc[a4'b]d4''^ + -4>^dab/4>d + ^ (j^ 9 c[aRb]'' 4> d " ^ ^[a^^bjc 
-3 f][a [06] (/-c - ^ (/-blc + ^ </'^-Rfc]c + ^ 4>^gb]cR - ^ Sbjc ^''^d] 



Unfortunately the terms introduced into the system by the scalar field, Tab and mabc, complicate the equations a lot 
due to the form of the energy momentum tensor for the conformally invariant scalar field. 
Furthermore there is an additional equation, 

f72i? + 6f7V"Vaf7-12(V"f))(Vaf7) = 0, (4) 

which must be satisfied at one point to be automatically fulfilled everywhere. 

If we fulfill equation (||) at one point the function i?(t, r) can be given freely. It determines the gauge freedom in the 
conformal factor to a certain extent. If [M^gab,^) is a solution of the unphysical equations so is {M,9'^ gab,dQ) 
for 6* > 0. In all the calculations, R has been set to 6, the value obtained if the compactification for Minkowski space 
given m ^ is used. The program code does allow to specify R{t, r). The conformal gauge freedom has been discussed 
m more detail in [||. 
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C. Simplifications by tiie splierical symmetry and tlie remaining gauge freedom 

The first order system is an underdetermined system. To make it complete the coordinate and frame gauge 
freedom must be fix. This will be done in this section. 

Furthermore the number of variables will be reduced significantly by making use of the spherical symmetry. 
A spacetime is said to be spherically symmetric if it possesses a isometry group / which contains a subgroup R which 
is isomorphic to the 3-dimensional group SO (3) and whose orbits are orientable submanifolds of dimension < 3. The 
orbits of R are either (fix) points or spheres The fix points form at most two timelike geodesies as can be seen 
by symmetry arguments. It is assumed that the initial hypersurface intersects with at least one line of fix points, 
i.e. there is a (regular) center on the initial hypersurface. 

But no attempt was made to use every simplification spherical symmetry provides for the following reason: It is 
well known that Einstein's equations can be reduced to two ordinary differential equations (constraints) plus 
matter equations in spherical symmetry. The circumstance that no equation with a time derivative for the geometry 
variables must be solved can be viewed as expressing the lack of gravitational radiation in spherical symmetry. I have 
not made any attempt to incorporate this kind of simplification into the system for the unphysical spacetime, since 
the calculations are supposed to be a playground for testing and learning about the advantages and disadvantages of 
the formalism for a later application to models with less symmetry. 

In the gauge used the system is symmetric hyperbolic and the coordinates cover the whole domain of dependence of 
the initial hypersurface. It is not known to the author how to construct such a coordinate system in less symmetric 
spacetimes. The very nice feature of semilinearity is certainly an artifact of spherical symmetry. 

1. Gaussian gauge 

In the analytical analysis of the initial value problem with scalar fields as source terms, a Gaussian gauge was 
used Q. This gauge gives a symmetric hyperbolic subsidiary system of equations. If the energy momentum tensor 
fulfills certain conditions, the coordinate system will break down because of the formation of caustics |]l6|. Whether 
the region of the unphysical spacetime representing the physical spacetime can be covered by the use of a Gaussian 
gauge in unphysical spacetime is not obvious since the conformal rescaling factor U, acts as a kind of artificial energy- 
momentum source in unphysical spacetime. By numerical calculations it was straightforward to show that even in 
unphysical spacetime a Gaussian gauge will lead to caustics , and is thus inappropriate to analyze the strong field 
regime. 

The system of equations is very similar to the one obtained in "double null" coordinates except that it is not semilinear. 
It will not be discussed any further. 

2. "Double null" gauge 

The construction of the coordinate system and the frame are described starting from an initial hypersurface. The 
remarks about the differentiability class of the objects assume that the construction is done with respect to a given 
C°° manifold (for a definition see section 1.1.1]). The investigation of the differentiability is necessary to ensure 
that a discontinuity of a value or even a singular value of a variable at the center is not caused by lacking smoothness 
of the coordinates. 

The gauge realized in this way turns out to be an obvious adherent of double null coordinates. 

Assume we are given a spacelike C°° hypersurface S with t = t^. This hypersurface is factored by the orbits of the 
group. By the geodesies running from the center in the different directions the angle coordinates (i?, Lp) with line 
element d'd^ + sin^ d dip on the unit sphere are defined on all orbits. By isotropy the geodesies are perpendicular to 
the orbits. The orbits are labeled by a monotonically increasing coordinate r, defined to be at the center. Auxiliary 
coordinates (u, v) are defined by u := (to — r)/2 and v := (to + 'r)/'2- Set (u, i3, ip) ~ const along the future directed, 
outgoing null lines, (f,??, (^) — const along the future directed, ingoing null line. On passing through the center set 
u — V. Due to the spherical symmetry there cannot be any null caustic, a break down of the coordinates is aligned 
with a spacetime singularity. By (w, v) a timelike coordinate 

t^u + v (5) 

and a spacelike coordinate 

r = V — u (6) 
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are defined everywhere in the domain of dependence of S. An orthonormal frame field (eo°, ej^", e2°, 63") is defined 
except in the center (polar coordinate singularity) by normalizing the orthogonal vector fields {dt"' , dr"' , d^'^ , d^"') . 
As shown in |^8| the t = const hypersurfaces are at least and hence sufficiently smooth. The proof proceeds in 
two steps, firstly it is shown that there is a C°° transformation to radar coordinates, secondly a theorem by U. ProfF 
| p9| , theorem 1.2.4] about radar coordinates completes the proof. 

In this coordinate system the following relations hold: The frame-coordinate matrix is diagonal, 

/eo"(i,r) 

ei}{t,r) 



and 

All Ricci rotation coefRcients except 

l-oi{t,r), 
7-22 (i,r), 

7-22 (i,?-), 



V 



e2'(t,r) 

e2^(t,r)/sin?? . 



7-33(^,0 =7-22 

7-33 (^''") = 7-22 



(7) 



(8) 



7-33 = 



cos ■& 
sin?? 



62 



(9) 



vanish. 

Scalars invariant under rotations are functions of (t, r) only, rotationally invariant vectors V"' are of the form 

/V^t,r) 
V^t, r) 


V 



(10) 



and symmetric covariant 2-tensors Sab, e.g. the Ricci tensor Rat, invariant under rotations, look like 



/ Soo{t,r) Sgiit,r) 

Soi{t,r) Snit,r) 

S22it,r) 

\ S22it,r) 



(11) 



This follows from the assumptions about the symmetry. 

All components of dabc"^ are either zero or proportional to dioi-- 



3. The resulting system 

Due to the complicated form of the matter terms the set of equations is very lengthy. The equations are given in 
appendix they are derived from system (^ and 

K^)-.(^)-^ 

The system (|^) has some remarkable features which should be mentioned. Firstly it is semi linear. Secondly every 
equation has one of the following forms: 

duU = buiU,V,T) 
d,V = bv{U,V,T) 
dtT = br{U,V,T), 

where U, V and T are variables propagating along u, v and t. lA stands for 78, R\ and (/>!, V for 78, RZ and 03, T 
for Co", e, 72, 7, i?2, fi, VLq, ili, (j), ^o, and dipi-. For every T there also exists a constraint, 

drT ^ ct{U,V,T). 

There are no constraints for the U and V. 
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4- The regularity conditions 



Polar coordinates cause regularity conditions in the center, in physical spa cetim e as well as in unphysical spacetime. 
For the variables in the system ( |B3| ) those conditions are given in appendix [B 3[ They express that locally the center 
behaves like Minkowski space, and indicate how a wave hitting the r = (inner) boundary is reflected there (passes 
through the center) . The regularity conditions at the inner boundary are the part of the code that would change if a 
spacetime with a throat were calculated. 

But in unphysical spacetime there are also necessary conditions for the regularity at — expressing an appropriate 
fall-off in physical spacetime and therefore reflecting asymptotical flatness. 

It is not known whether those conditions are also sufficient for regularity on in the general case. The answer to 
this requires the investigation of the constraints on a given hypersurface. For the case that a neighbourhood of J is 
free of matter and for a special choice of a hypersurface (for technical reasons) it has been shown that the necessary 
conditions are also sufhcient |^^. In the calculations presented here those conditions are fulfilled. 
There are as many regularity conditions (at the center and at J^) as there are variables T. 



III. THE NUMERICAL METHOD 



In this section the numerical methods used will be described shortly and the reasons given for choosing those 
methods. 



A. The initial value solver 



There are 12 constraints and 12 regularity conditions to be solved. The free functions in the constraints are the 
six variables U and V. But giving those makes the interpretation of a parameter study difhcult. U and V represent 
higher derivatives of the primary quantities, the metric gab^ the rescaling factor 17 fixing the relation between physical 
and unphysical spacetime and the scalar field (p. Furthermore the regularity conditions at are easier to handle if it 
is known where f2 vanishes. A straightforward way to realize this is to give SI on the initial slice. 
In the code e]_^, 7-11, 7-22, SI, 0, and 0o are given as free functions (because of regularity they all must be even at 
the center), ei^ = eg" determines the 500 respectively the gn component of the metric, 7-11 and 7-22 the extrinsic 
curvature of the initial slice in unphysical spacetime. The constraints which contain derivatives of these quantities 
are interpreted as algebraic conditions for the quantities on the right h and side. To ensure the invertibility of the 
resulting system, including the points where S7 = 0, instead of equation (B29) its derivative is used. A constraint for 

R\ + Ri results which becomes an algebraic condition on SI = if Rl + R'i is at least (condition B38k ). This is a 
necessary condition for a regular S . 

The differential algebraic system has boundary conditions in the center and at = 0. As there are as many boundary 
conditions as variables to solve for all the initial value freedom is coded in the free functions ei^, 7-11, 7-22? </>, 
and 00 . The system is solved with a relaxation scheme combined with a Newton-Raphson solver derived with minor 
modifications from the code given in |2l]. Firstly the system is solved between r = and S7 = 0. If desired the same 
scheme can be used to extend the initial hypersurface beyond the intersection with . For this integration boundary 
conditions are given at S7 = 0, namely the values obtained at S7 = when solving the constraints "inside" . To avoid 
a dependency of the values at J on the treatment of the outer grid boundary the initial hypersurface has always be 
slightly extended beyond J. 

To improve accuracy the system is solved with different grid sizes and the results are Richardson (or Bulirsch-Stoer) 
pl| extrapolated to vanishing grid size. In the test cases, where an exact solution is known, the accuracy was limited 
by the rounding error of the inversion of the matrix in the Newton-Raphson part of the code. Typically the numerical 
and the exact solution differed in the last two digits for calculations with 8 byte floating point numbers. 



B. The time integrator 



Constrainted evolution schemes are known to give in general more accurate results than free evolution schemes. But 
there are disadvantages of constrainted evolution schemes. In a constrainted evolution scheme the values at a certain 
grid point depend on the values at many other grid points on the same time slice. Thus if one grid point becomes 
singular the property of being singular is spread over the time slice into regions which are not causally connected with 
the singular point. 
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In most approaches used in numerical relativity the singularity is avoided by an appropriate choice of the coordinate 
system, actually the necessity of singularity avoidance has become a dogma of numerical relativity A quite 

often used example for a singularity avoiding coordinate system is obtained by the radial gauge in spherically symmetric 
spacetimes ||5|,p4[]. This coordinate system cannot penetrate apparent horizons which are supposed to wrap around 
singularities by the cosmic censorship hypothesis. The very interesting region of spacetime, where the gravitational 
field has become strong enough to prevent light rays from expanding, is excluded from being calculated. 
The double nulllike coordinate system used does not avoid singularities, the coordinate lines can only end at a 
singularity. In this subsection it will be described how a program crash is avoided if singular points occur and how 
the calculation is continued in the region of spacetime which are outside the future domain of dependence of singular 
points. In the chosen approach the singular boundary of spacetime is represented by grid points. 

1. Inside the grid 

To avoid the spreading of the singular property out of the domain of influence of the continuous equations it is 
necessary to run the scheme, at least in the very neighbourhood of a singularity, with a Courant factor of 1 . Different 
schemes have been tested |Q . The only second order schemes which could be run at a Courant factor of 1 were the 
second order schemes of the class |^ with /3 = 1/2. A second order Lax-Wendroff scheme as follows has been 
chosen: For the equation dtfi = Xidrfi + b{f), with A/ — const and a vector of functions /, the predictor step is 

f^iji = \ (fi-. + + ^ (^^ - fi-^ + T K ^ ^-^ 

Then the corrector step is 

If Ai were a function of /, as it is the case for other coordinate choices, it would be discretized as b. 
The second order Lax-Wendroff scheme is used to evolve gridpoints not depending on the values at the boundary. 

2. At the boundaries 

The treatment at the outer boundary is irrelevant as long as the scheme remains stable for the following reason: 
The initial value surface is extended beyond the intersection with J^. Running with a Courant factor of 1 the points 
inside and on J^, representing the physical spacetime, do not even depend on the values at the outer boundary which 
has no intersection with J^. Even if the treatment of the outer boundary causes an instability it will not influence 
the physics, i.e. the values on M and J^. In most runs with Coraunt factor 1 the grid points depending on the outer 
boundary have not even been calculated. For a run with a Courant factor < 1 the influence of the outer boundary 
treatment on the physical spacetime must converge to with the same order the scheme converges. Thus as long as 
the treatment on the outer boundary is numerically stable the values at M and are in the limit of vanishing grid 
size independent of the outer boundary treatment at any physical time. 

Finding a stable treatment of the center was a very difficult problem. So far I could not find a treatment which 
remains stable if a fourth order scheme is built from the second order Lax-Wendroff by Richardson extrapolation. 
Especially it was not possible to extend the grid from the gridpoints at r = and r = Ar to r = —Ar and run the 
Lax-Wendroff scheme on the extended grid. 

The solution of the stability problems at the inner boundary was to replace the values near the center after every 
time step by the values obtained by other methods. 

In the method used for my Ph.D. thesis the constraints were integrated from grid point number 2 towards the center 
(the constraints together with the regularity conditions determine all the variables at the center). The values at grid 
point 1 have been obtained by interpolation. The solution did not look smooth on the innermost gridpoints for coarse 
grids in the very strong field regime. For very large fields (values of A beyond 1.10) this method even became unstable. 
Therefore another method has been developed. 

In the calculations for this paper a kind of polynomial extrapolation with dissipation has been used to replace the 
innermost gridpoints. The idea is the following: Use the values at gridpoints 2 and higher and the regularity conditions 
to extrapolate towards the center by polynomial fitting and get solution I. Do the polynomial fitting again starting 
at grid point 3 to get solution II. Solution I and II can now be added in such a way that the simplest grid mode 
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with values 1,-1,1,... is eliminated. This adding of dissipation is necessary to ensure stability. I call this method 
polynomial extrapolation with dissipation. 



3. The singularity catcher 



Since the coordinate system does not avoid singularities some variables may become singular. Depending on the 
default setting of the compiler this causes a crash of the program by a floating point exception (on UNIX systems the 
signal SIGFPE is sent). The programming language C allows to specify what to do in case of a certain signal. In the 
program used the action on SIGFPE is to flag the corresponding grid point as singular and to continue the calculation 
on the rest of the grid. In addition to the reception of a SIGFPE signal a grid point is flagged as singular if either 
the principal part, i.e. eo° or 1 — 1 17^ 0^, of the equation changes sign, i.e. the system of equations becomes singular 
in the sense of 0| , or if the evaluation of the values according to the scheme would involve points already flagged 
as singular (hence called influence singularity). According to the latter reason every point whose values depend on 
values at singular points must be flagged singular. 

A necessary condition for the stability of a scheme solving symmetric hyperbolic systems is the Courant-Friedrichs- 
Lewy condition. The domain of dependency of the discretized equations must be a superset of the domain of depen- 
dency of the continuous equations, the Courant factor must be < 1. A Courant factor of 1 means that both domains 
of dependency coincide. On the other hand, every grid point depending on singular gridpoints is singular. Thus all 
points in the future null cone i of a singular point must be flagged singular. If the Courant factor is < 1 points outside 
L also. For a scheme with a Courant factor of exactly 1 at least in the very vicinity of the singularity it is possible 
to distinguish timelike/nulllike singularities from spacelike singularities. If a line of gridpoints is numerically singular 
only because of a influence singularity it is a strong hint for a timelike singularity behind. Pure influence singularities 
appeared only due to instabilities during the tests of different treatments of the center. This is in agreement with 
cosmic censorship. 

The dependence of the position of a singular line on the grid size has been tested intensely. In all calculated cases 
presented the position changed only by a negligible amount corresponding to a few grid points for calculations with 
about 10000 spatial points. The dependence of the first appearance of a singularity on the Courant factor has also 
been tested. It only changed by a few grid points when using Courant factors significantly smaller than 1. 



On the exact solutions given below a number of tests have been performed with grids varying from 300 to 10000 
spatial grid points and different Courant factors. For most plots a coarser grid has been used, typically with 100 x 100 
grid points. 

On the plotting grid the scheme shows the expected convergence behaviour. For finer grids (2000 spatial points and 

more) the error is dominated by second order terms. For coarser grids fourth order terms (oscillations corresponding 

to high Fourier modes) dominate the error near regions with steep gradients. 

The error predictions by Richardson's extrapolation are in good agreement with the actual error. 

The violation of the constraints converges at least with second order. Convergence is partly dominated by the error 

in the discretization of the derivative for an evaluation of the constraints. The violation of the constraints turned out 

to be an inappropriate measure for the quality of the solution. 



A scalar field on a flat background is obtained by setting k = in the evolution equations. Since the physical 
energy-momentum tensor Tab is tracefree the Ricci scalar R in physical spacetime vanishes. (f> is determined by 
□ — 0. The corresponding solution in unphysical spacetime is (p — (p/il. On M the following rescaling has been 
used: 



IV. CHECKING AGAINST EXACT SOLUTIONS AND ACCURACY ESTIMATES 



A. Comparing with exact solutions 



1. Scalar field on flat background 



±2 



(12) 
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with an obvious choice of the sign, depending on the position relative to J7and i+. 



^-i/(i-.^(.-.of) \r-ro\<l/a 



/M= ' ' for l'---o|^,V- , (13) 

1^ otherwise 

with constants a and ro- 

This solution is useful for testing the matter part and the geometry part of the code separately since k = decouples 
their evolution equations. But as most of the geometry variables are just constants this is only a first step as test for 
the geometry part . 



2. By conformal gauge freedom deformed background 

For any positive function '&{t,r) the spacetimes {M,'d'^gab,4>/'&) ^^J^d (Ad,gab,4>) correspond to the same physical 
spacetime if {M, gab, 4") is a solution of the unphysical equations. This gauge freedom can be used to obtain a 
unphysical representation of Minkowski spacetime where the geometry variables evolve in time and are not constant 
(except the Weyl tensor component dipi-, which is zero for every representation of Minkowski space). For 6 an ingoing 
wave pulse with the initial form ( p^ ) plus a constant c, such that 9 > everywhere, has been chosen, (f) has been set 
to 0. 

This solution has been used to make the geometry variables vary in time and space and is thus a better test for the 
geometry part of the code than the test in the previous subsection. 



3. (j) + const solution 

In physical spacetime 4> = const is a solution for Minkowski spacetime. If we use const (p = is still a 
solution on M in unphysical spacetime. But this solution is no longer regular on J' thus the calculation can be done 
in the interior of M only. For the tests that poses no restriction, since in the unphysical spacetime there is nothing 
special in the evolution equations on J^. 

Since in this solution k ^ and <p =^ G the coupling between the geometry and the matter part can be tested. 



B. Error estimates and problems near the critical parameter 



For the nonlinear regime two ways have been used to check the accuracy of the solution. 
The first is the Richardson extrapolation to vanishing grid-size with the resulting error estimates. The second uses 
the fact that the position of J is known and thus it is possible to compare the calculated value O there with the exact 
value 0. The error estimate for 17 by the Richardson extrapolation was in all cases in agreement with that deviation. 
On an approach to the critical parameter, the smallest parameter where a singularity appears, there is a dramatic loss 
of accuracy although on the coarse plotting grid the scheme still converges quadratically. Viewing the fine grid used 
for the calculations one sees that the loss of accuracy is largest around gridpoints 5-10, and this maxima of the error 
decreases linear with the grid size, meaning that there is only quadratic convergence in a norm, but not pointwise. 
For models not to close to the critical value the accuracy problem described can be cured by using more gridpoints. 
But with about 50000 spatial gridpoints that brute force method reaches its limitation. Since the scheme converges 
quadratically in a norm only, building a higher order scheme by Richardson extrapolation did not work. 
A detailed inspection of the reasons hints, that this is partly a problem of numerical regularization as described 
in j26j (pages 13-15). When the calculations for this paper were done there was no way known to the author to 
numerically regularize the Lax-Wendroff scheme. 



V. CALCULATIONS 

In the calculations presented the free functions defining the geometry of the initial slice are given as follows: 

ei^ = 1 (14) 
7-11 = (15) 
7-22 = (16) 
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±2 

, (17) 
^'(l + tan2 *2±r)(l + tan2 ^a^^) 

The value of the time coordinate on the initial slice tg is 7r/2. The form of is an often used choice for the rescaling 
of Minkowski space. 

For the scalar field the initial value is given by 

otherwise 
with (7 = (S/tt) and ro = 7r/4. 

This function 0(r) is the uniquely defined C^-function with compact support [ro — 1/ct, + l/tr], being a polynomial 
of degree 8 on the support and with maximal value A a,i tq. A is the parameter to be varied. (j){r) should be at least 
in order to have all variables in the system at least avoiding the known problems of a Lax-Wendroff scheme at 
discontinuities. From the numerical viewpoint this form of the pulse is much better than a C°° partition of the one 
like ([T^), since the variation of (j^ir) and its spatial derivatives is better distributed over the support. 
(pQ is chosen in such a way that the pulse would be purely ingoing on a flat background, which is the Einstein cylinder 
in unphysical space, i. e. 

, cosr , 
^o = (t>i+(t>- — (19) 
~ ~ sm r 

for r — ro G] — 1/cr, l/cr[ and otherwise. 

A pulse with compact support provides several advantages when checking the accuracy and interpreting the numerical 
solution since part of the qualitative behaviour of the structure in the large is known from analytic considerations. 
Figure ^ shows the qualitative behaviour in the large for a regular solution. Region I is the compact support of the 




ingoing pulse in the linear model (k = 0), passing through the center and then crossing the thick line. Every 
deviation from a Schwarzschild solution in region III is caused by a pure back scattering effect. The deviations from 
flat space in region II are also an back scattering effect — the field gets caught near the center for some time. As 
the characteristics have slope 1, region IVa is a portion of flat space, region IVb is a portion of Schwarzschild. For 
stronger data the solution will become singular — before reaching the tip a singularity will appear. To investigate 
the structure of this singularity is one of the goals of this work. In a regular spacetime timelike infinity, , lies at the 
point (tt, 0). 
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A. A parameter study 



With increasing parameter A more and more matter is put onto the initial shce. At some critical parameter A* the 
spacetime is expected to become singular. In this subsection the change of the global structure will be discussed and 
a conjecture for the critical case presented. 

When interpreting the behaviour with varying parameter two things have to be kept in mind. 

Firstly the characteristics, null geodesies, are a common structure of all models. As the coordinates are adapted to 
this structure it is straightforward and well defined to compare the models. 

Secondly the area distance f := ^^ggg to the center in physical spacetime, expressed in unphysical variables, is 

r 

Even on the initial slice e is determined as a solution of differential equations, the constraints. Depending on A the 
area distance of the inner and outer edge of the shell of compact support of (p vary (in this parameter study they 
monotonically grow with A). In principle it is possible — and for a better comparison with M. Choptuik's work 
desirable — to give e as a free function, but parts of the initial value solver code would have to be rewritten. But 
in doing so one can no longer give either Q or components of the extrinsic curvature as free functions. The first case 
causes problems since the location of 5 H on the grid is no longer known, in the second case it is unknown whether 
the necessary conditions for regularity on are also sufficient. 

1. The conformal factor 

Figure shows a contour plot of the rcscaling factor fl for modestly strong initial data corresponding to a parameter 



n(t,r), A=0.25 
3_5 K ■ 




0.0 0.5 1.0 1.5 

r— axis 



FIG. 4. n{t,r) for A = 0.25 

A — 0.25. That is well below the critical value A*, which is somewhere in the interval ]0.48,0.49[. Apart from 
some minor deformations in the contour lines of fl, the figure looks like Minkowski, although nonlinear effects like 
backscattering are already significant. The 17 = contour coincides with the analytic location of J' represented by 
the thick dashed line. The values depending on the outer boundary of the grid have not even been calculated. That 
is why there are no contours in the upper triangle. 

In figure ^ the spacetime has become singular. The thick line is certainly a singular boundary of the future domain 
of dependence of the initial data since here and in all the following cases scalar invariants become singular. Near the 
center the singularity is spacelike. Further outside its slope is numerically indistinguishable from a null line. The edgy 
look is an effect of the plotting program and the coarse plotting grid, which will be seen later when the singularities 
will be examined more thorough — a finer grid near the singularity showing more details will be used there. 
With increasing amplitude of the initial scalar field (figure ^) the spacelike part of the singularity has grown. 
The change of the shape of the singularity on approaching the critical parameter from above suggest the following 
conjecture: In the critical case the singularity is a null line. 

For further increasing A the outer edge of the singularity becomes spacelike. This part grows inwards and finally 
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Q(t,r), A=0.49 



J. '3 



3.0 



2.5 



2.0 
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r— axis 



FIG. 5. 0(4, r) for A = 0.49 



[]('t,r), A=0.55 




FIG. 6. n{t, r) for A = 0.55 



n(t,r), A=1 .5 
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FIG. 7. Q{t, r) for A = 1.5 
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we get a picture like figure |^. Note the contour lines of fl are focused at the intersection of J7_with the singularity. 
Outside the scalar il goes to — oo, inside to oo on approach to the singularity. In subsection V A2 the structure of 
it will be examined. 

The critical case can also be approached from the subcritical side. In figure || the time dependence of Q in the center 




is plotted for the three subcritical values A = 0.25, 0.45, 0.48 and the supercritical model A = 0.49. 
Approaching the critical case O(i,0) seems to develop a zero point before «+. For spacetimes above the critical case 
the value of ^l{t, 0) at the singularity is greater than but again approaching on approach to the critical case. 
Outside the center and away from J', the scalar ^l does not seem to go to even in the critical case. Although 
the values of f2 are influenced by gauge, the statement that f2 vanishes is not. The author cannot interpret this is 
exceptional conformal structure for the critical case yet. Nevertheless the calculations show the critical model might 
be recognized by this behaviour of fl. 



2. Are the singularities found of physical nature? 



The figures of the preceding subsection show that there are singularities for large values of A. In subsection |111 B 3 
it has been described when and how the program flags a point as singular. That does not necessarily mean that those 
points represent a singularity in physical spacetime. In this subsection arguments will be given to decide that issue. 
Furthermore the quantities used classify the singularity. 

In their singularity theorems R. Penrose and S. Hawking have shown that a singularity is unavoidable if there are 
trapped surfaces and the energy-momentum tensor fulfills certain energy conditions. The conformally invariant scalar 
field does not fulfill those energy condition. Nevertheless I regard the appearance of trapped surfaces as a strong hint 
for singularities. With the transformation to the solution of a massless Klein-Gordon field it is easy to check 
for the existence of an apparent horizon in a massless Klein-Gordon field model which satisfies the relevant energy 
conditions. 

The null expansion of in- and outgoing null directions in spherical symmetry can be written as 

^~o„M„ =^ (s„,."V,r), (20) 

where is the outgoing null vector and 

e." = (eV-eV) 
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is the ingoing null vector. Written in unphysical quantities: 

^o„t = ^ (7-22 -l)~^o-^i 

= fl (7^22 + ^)-^a + ^i (21) 

There is the freedom of null boost, thus 9^^^ and 9^^ are not gauge invariant, but their product is. On the first view 
one might expect that 9^^^ vanish for r —^ 00. The gauge for ^ chosen makes 9^^^ not vanishing on j/"*", 9^^ not 
vanishing on J^^ . That is a pure gauge effect, by an appropriate null boost 9^^^ and 9.^^ can be made vanishing on 
J. 

Another criteria for or against a "real" singularity is the behaviour of curvature invariants of physical spacetime, e.g. 

Cabcd.C"''"''^ = fl^ dabcd d"''"''^ — 12^1^ (c^l01~) • 

Figure |^ shows the spacetime region near the singularity of a A = 0.55 model. The thick, dashed line is a part of 

6„„„ A = 0.55 

2.80 
2.70 
m 2.60 

D 

^ 2.50 
2.40 
2.30 

0.00 0.10 0.20 0.30 0.40 0.50 
r-Qxis 

FIG. 9. Bout for A = 0.55 



9„„„ A = 1.5 



2.6 
2.4 

O 

I 2.2 
2.0 
1 .8 

0.0 0.2 0.4 0.6 0.8 

r-axis 

FIG. 10. Sout for A = 1.5 

the thick line the singularity S. The singularity is covered by a region of trapped surfaces, the thin line shows the 
apparent horizon. In the corresponding initial value problem for a massless Klein-Gordon field the same qualitative 
picture arises. 

The proper time of a Bondi observer (see equation (p5[)) goes to infinity on approach to the intersection point V of 
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the S with This is a strong hint that P is a singular timeHke infinity, a singular i+. Thus the last null line from 
the center which reaches J'^ at i"*" is the event horizon (dotted line). 

These are already strong hints that 5 is a singularity of the physical spacetime. It is further strengthened by the fact 

that the curvature scalar I Cabcd C"''^'^ j blows up at least on approach to the inner (spacelike part) of S. There are 

also hints that the area radius goes to on S, but to decide that issue more numerical accuracy near the singularity 
would be needed. 

All these considerations "prove" that 5 is a real singularity of physical spacetime. 

The situation is quite different in figure Along S the scalar il goes to +oo inside J""*", and to — cx) outside 

S is spacelike everywhere. (Cabcd C"''"^'') does not show any sign of blow up near S, there is no region of trapped 
surfaces wrapped around S. Due to the discontinuity of near V the formula for the proper time of an observer at 
is not applicable. 

Since only quantities of the unphysical spacetime become singular on 5 I conclude that 5 is a conformal singularity 
caused by the gauge chosen by specifying R{t, r) — 6. It is very similar to a coordinate singularity. Only minor, 
unsuccessful effort has been spent to try to avoid the conformal singularity by using different choices of R{t,r). 
There is another very interesting point in the models with A > 1.2. On the initial slice there are regions where 
both 0.^ and 9^^^ are positive ( "antitrapped surfaces"). A calculation back in time shows a singularity in the past, 
covered by an antitrapped region, i.e. a spacetime with a white hole. Although the model A = 1.2 shows a conformal 
singularity in the future the existence of an apparent horizon suggests that the physical singularity is hidden behind 
the conformal singularity. It is natural to assume that beyond A* there is always a physical singularity which is 
sometimes hidden by a conformal singularity. 



3. The mass scaling relation 

In his study of scalar fields M. Choptuik found the critical mass scaling behaviour msH — a{A~ A*)'' with 
independent of A and an exponent 7 of approximately 0.37. 
The Hawking mass 

m(i,r) = ^ {l + P'e^^J 

can be written as 



1 /r\^ 



0' ( (7-22)' - (J) ) + 20 (7^22 '/'o + J0i) + ((0o)' - (0i)') 



-l-K(f>^n. (22) 

8 e 



The use of equation (B36) on is essential to write the Hawking mass in a form not containing terms with 
factors. On the Hawking mass is identical with the Bondi mass. For scalar field data with compact support on the 
initial slice in M the Bondi mass on the initial slice coincides with the ADM mass. 

Figure |ll| shows an logarithmic plot of the black hole mass obtained for the parameter study performed here (the 
stars 7k- mark calculated models). The black hole mass is determined as the Bondi mass "at i+". 
The value of 0.4866 for the critical parameter A* has been found by optimizing the straight line look in figure [ll|. 
The critical exponent is 0.37, but due to the uncertainty in A* the error is certainly ±0.05. Hence the results are 
compatible with the results of others, but I would not claim more. 

Although the mass caught by the black hole increases the percentage of the caught mass has a significant maximum 
at A w 0.75. For large A almost all the mass escapes to Actually A >w 1.0 the given value is only an upper 
bound since is hidden behind the conformal singularity and thus the Bondi mass at the intersection of JJ^ with the 
conformal singularity may further decrease. The figures for the black hole mass of the corresponding Klein-Gordon 
models show the same behaviour. 

For models with large initial scalar field amplitudes A there are number of significant changes in the structure of 
spacetime. Although most of the scalar field is still contained in region I the shell with most of the Hawking mass in 
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it expands and crosses J in region III. Mass and scalar field "decouple" , propagation of mass is almost completely 
determined by nonlinear terms. 



B. Christodoulou's theorem 



D. Christodoulou investigated the general relativistic massless Klein-Gordon field in spherically symmetric space- 
times in a number of papers, for a list of references see p7|] . He found sufficient conditions for the appearance of a 
spacetime singularity. The theorems proven there contain the following statement: 

Theorem 1 Consider on the initial future null geodesic cone Cq an annular region bounded by the two spheres 

Sifl and §2,0 with §2,0 *^ the exterior o/§i,o and areal radii ri_o and ^2,0- Let 

So-.^'r^-i e(0,i), 



Vo 



2 (m2,o - mi,o) 



r2,o 



E{y) ■■= 



log 



2y 



5 — J/ 



A sufficient condition for a non-timelike singular boundary in the future of Cq is 

Vo > E{6o). 

The model with a conformal scalar field calculated here is (almost) equivalent to the massless Klein-Gordon field 
[^. If (M,5^j, 4>) is a solution for the general relativistic conformal scalar field model {M,gab, (j)) is a solution for the 
general relativistic massless Klein-Gordon field with 



4> — \/6arctanh— 

, 1 ~2, 

9ab = (1 - -0 )gab- 

The area radius f and the Hawking mass in are given by 

m = m + — I 0^^^ [lo) + 6^^ (uj) H (oj) e^^ (uj) 



(23) 
(24) 



1 — if22(/)2. f is the area radius, m the Hawking mass, 6^^^. the null expansion of the out- and ingoing null 

directions with null vectors e^ and in the conformal scalar field model, fh can be written in a form which does 
not, even implicitly, contain terms proportional to fl^^. For the purpose of this subsection it is not necessary to write 
m in a form obviously regular on J^. 

In figure 13 r/o — E{So) is plotted for three outgoing light cones versus the area radius f. For a parameter oi A = 0.75, 
a model far beyond the critical value with a distinct singularity in physical spacetime, the solution has been printed 
out for 100 time slices with 100 grid points. On outgoing null cones 770 — E{6o) has been evaluated. The sufficient 
condition is satisfied if 770 — E{6a) > 0. The gridpoints are marked with squares, diamonds, and triangles. The first 
light cone (squares) lies outside the event horizon, the criteria is nowhere fulfilled. The second lightcone (diamonds) 
approximately coincides with the event horizon. The criteria is not fulfilled either. Only if the apparent horizon is 
crossed (triangles) 770 ~ E{do) becomes larger than 0. Although many outgoing lightcones in various singular models 
have been checked the criteria was only fulfilled when the lightcone crossed the apparent horizon. Thus I conclude 
that the criteria is not very sharp. This is in agreement with the results in ES]. 
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FIG. 13. Christodoulou's criteria rio — E{So) versus area radius r 
for various light cones 



C. Extraction of Radiation 

In addition to the already presented determination of the Bondi mass I am going to demonstrate the simplicity of 
radiation extraction on two more examples in this subsection. In the first case the effects of the nonlinearities in the 
equations are analyzed — the purely backscattered radiation in region III and the decay of the radiation in region 
IV of a regular spacetime model. In the second case values of the Bondi mass in a singular model are compared with 
values of the mass read off at finite radii. 



1. Effects of non-linearity 

Figure |lj compares (j) on J', i.e. the coefficient of the 1/r term of the scalar field in physical spacetime, for a 



Scalar field on S 
0.2 I 




FIG. 14. (j) for the whole time evolution 
model with an initial amplitude oi A = 0.40 for the linear model {n — 0, dotted line) and the nonlinear model (k = 1, 
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continuous line). In the model with gravitation a significant amount of radiation, reflected outward by backscattering, 
has already crossed J before the linear signal reaches J (region III) . The main signal is stronger in the nonlinear case 
and there is some scalar field left in region II. 

In the figure |l^ the Bondi mass is shown in the pure backscatter region III. Later times correspond to the backseat- 
Change of Bondimass in recior II 
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FIG. 15. rriBondi in region III 
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FIG. 16. in region II 



tering of a matter shell which has decreased in size during the infall. 

The scalar field </> is displayed in a double logarithmic plot of the scalar field versus the proper time r of an observer 
at J . This observer is obtained by taking a world line at fixed area radius f (and angle) and taking the limit f — > oo. 
In the unphysical variables one gets 



To = 



(25) 



where r j (t) is the coordinate r in unphysical spacetime of at unphysical time t and 



(eo°) 



o^2 



R-6R2 
12 



1 



+ ^ (eo" + 7) - 

rjj-^ Co 



(72)^ 
1 



'J 



(eo" + 7)' 



(26) 
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The dots represent grid points of the calculation with 5000 grid points. Since on approach to i+ the center is 
approached and the values become inaccurate, only those points have been included in the plot where the results of a 
Richardson extrapolation with 5000 and 10000 grid points are visually indistinguishable from those of an extrapolation 
with 10000 and 20000 grid points. In this plot the last 10 (of 5000) points on J before i+ are missing. 
The slope in the asymptotical region is approximately 0.34. Thus for very late times ~ (r — to)^ '^*. 
Unfortunately the Bondi mass falls off so fast that rounding errors of the graphic program, with which the calculation 
of the mass is done and which has only single precision, do not allow to enter the asymptotical region with the constant 
slope for 0. 

But in spite of the accuracy problem near the center, figure |l^ illustrates an advantage of the conformal method: The 
time scale of the collapse of the initial matter shell is of order 1, the fall off of the decay could be investigated for 
almost 10^ times the dynamical time scale. 



2. Radiative quantities at finite physical distance and at J 



A rough estimate of the errors which are to be expected by reading of radiative quantities at finite radius has been 
made in the following way: Given a lower bound for the area radius, the grid on a time slice has been searched for the 
first grid point P with larger area radius. At that point P the value for the Hawking mass ?TT.Haw,fin and the "finite" 
Bondi mass, 

1 / \ 

mBondi3n(-P) :== 2 (gj ^101" ' 

have been calculated. Those values have been compared with the Bondi mass on the intersection point of the future 
directed Hght cone of P and . The model used corresponds to the parameter A ~ 0.55 and develops strong 
gravitational fields, the spacetime has a singularity which catches approximately 85% of the ADM mass of 0.22. The 
outer boundary of the compact support of the scalar field lies at an area radius f = 2.7, P has been chosen two, five 
and ten times that value, expressed in units of the ADM mass, P is at 10, 25 or 50. 

The maximal relative error, i. e. error over ADM mass, is approximately 2%, 1% and 0.1%. But the relative error 
in the backscattering region III (error/mass loss in region III) is by a factor of 10 larger. Those errors are the errors 
made by approximating null infinity by grid which ends at a finite physical distance. There is almost no difference 
between mHaw,fin and mBondi,fin- 



VI. SUMMARY 



In the work reported in this paper a mathematical formalism, designed to investigate the structure of a spacetime in 
the large by conformal techniques, has been introduced into numerics. This paper demonstrates that global properties 
of a spacetime are numerically calculable and must not necessarily be estimated from a numerical simulation by 
heuristic arguments. 

The set of equations can be solved either by a spacelike or a characteristic initial value problem. The form of the 
equations differs from Einstein's equations. The technical tricks developed for the numerical solution of Einstein's 
equations cannot be used unchanged. The accuracy problem as described in section IV B shows that there is still the 
potential for further progress. Recently, too late to be included in this paper, a new way of discretizing was found, 
which propagates the regularity condition at the center without switching schemes near the center, as done here. The 
accuracy has improved by a factor of 10 for models near the critical parameter. 

In this formalism important properties of spacetimes are automatically fulfilled, the ADM mass is always conserved, 
the energy radiated away equals the change in the Bondi mass, and there is no spurious ingoing wave caused by an 
unphysical reflection at the outer boundary. Furthermore the location of event horizons is straightforward and the 
calculation yields "Penrose" diagrams allowing to read off the causal structure of singularities. 

The obvious advantages might not be worth the effort for the investigation of astrophysical questions where most of the 
uncertainty originates in an approximate equation of state. But if one is interested in the use of numerical relativity 
as a substitute for the lacking experimental approach to the open problems of mathematical relativity, the conformal 
techniques provide a very promising approach. In this circumstance it should be mentioned that H. Friedrich has 
developed the conformal techniques further |2^. There the rescaling factor is no longer a variable of the system 
fixed in a complicated way by the gauge source function R{t,r). The gauge freedom in the rescaling is fixed directly 
by specifying Q{t,r), and the problem with spacelike infinity may be solved ^o| . 
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APPENDIX A: NOTATION 

The signature of the Lorentzian metric gab is (— , +, +, +). 
Whenever possible I use abstract indices as described in pl^, chapter 2]. Small Latin letters denote abstract indices, 
underlined small Latin letters are frame indices. For the components of a tensor with respect to coordinates small 
Greek letters are used. The frame (^§77) is constructed from the coordinates a;'', e^" denotes an arbitrary frame. In 
this notation Va is a covector, a scalar, namely Va ei°. 

v{f ) is defined to be the action of the vector v"" on the function /, i.e. for every covariant derivative V a'- t{f ) = t°- V af ■ 
The transformation between abstract, coordinate, and frame indices is done by contracting with e^" and e^^. All indices 
may be raised and lowered with the metric qab and the inverse g^^ . g^'^ gcB — 5^b, A and B are arbitrary indices, 
e.g. Cia = gab &l and e^a = g^- eta- 

For a frame Ci" and a covariant derivative Va the Ricci rotation coefficients are defined as 
From this definition follows 

et''eh{^at')^e4tl)+ji^A 

With respect to a coordinate frame e^" = (gfji)" the components are the Christoffel symbols F"^^!^. 
The torsion T°'bc is defined by 

VaVb/- VfcVa/ = -T'abVJ, 

the Riemann tensor Rabc^ by 

VaVfcWc — Vb^a<-^c — Rabc^ — T^ab'^ d^^c- 

Contraction gives the Ricci tensor, 

Rab — R-acb i 

and the Ricci scalar 

R — Rabg'^'^- 

The Einstein tensor is given by 

Gab = Rab — 2 Rgab- 

The speed of light c is set to 1. The gravitational constant k in Gab = uTab has been set to 1 or in the calculations, 
corresponding to the full nonlinear theory or a scalar field on a flat background respectively. 

APPENDIX B: THE SYSTEM OF EQUATIONS IN DOUBLE NULL COORDINATES 

1. The time evolution equations 

From the system (|^) the following set of equations can be derived, where the abbreviations 
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and the substitutions 



71 = 


7-01 + 7-11 


72 = 


7^22 




73 = 


7^11- 


7% 


m = 


(-R01 - 


V i?ii)/2 


R2 = 


Roo — 


Rii 


m = 


{Rii- 


- i?oi)/2 




(001 + 


0ii)/2 




(011- 


0oi)/2 






R 






6' 



and 



-Oi?2--f23 (Too-Tii)-72f2o- -fii, 
4 4 — — - J. - 



<^2 = ,^00 - 011 ^ -272 00 - 2 J .^1 - ^ 



which follow from spherical symmetry, have been used. 

(5o+ai)73 = -{^^^l)+dwr^+\{^-R2) 

{dQ + di)m = (dioi-fio + f^m + %3)/(l-K (f7 0/2)' 



+ ( - 7I i?2 + 2 ^) + 72 (iTl - :R3) + ^ (5ii?) - ^ (aoi?) ) 



with 



7 ( (1 - K (17 ,/./2)^) (M + - i m)/r 
-f 17(01 /r) (17o0 + 0of^)) 



for r ^ 



§7^(^0 + 00^^)5^01 forr = 



(^0 + oi) 03 = 72 {(1)1 - (^3) - 2 71 03 



+ ^ (^-(m-i?3)-rfioiaO+li?2-ii? + 2727l 

+ ^(m + ^3+^i?) + (5ii?)+l7lii-l (aoi?))+5; 

7 ( ((0iA) 7 + 72 00 + ^ 0i? + 01 + 03) / 
+ 7I (0i/r)) 



for r 7^ 



771 9r0]_ for r = 



(ao - ai) 71 = - (73 71) + rfioi^Ji + ^ (f - -^2) 

(5g-ai)m = (dioi^ no + ^^m + 5^,) /(l-K (170/2)2) 
+ (-73(^^ + 2M)-72(m-:R3)-^ (^li?) - ^ (aoi?) ) 



24 



- f 0(<^i/r) (Oo'A + '^of^)) 



for r ^ 



for r = 



(ao-9i)(/)l = - 72 (^-^)- 2^73 

+ ^ (^m-:R3-dioi^f2+^^-^i? + 27872^ 



^ ( -(i?l + i?3)- ^i? 



~i M + ^^3i?--i(aoi?))+5; 



01 



>1 



7 ( - ((-^1 A) 7 + 72 (Ao + ^ <^ i? + + <^3) A 
+(<^i/r)73) 



for r 7^ 
for r = 



ieo°(73 + 7l) 



9oe ^ 



672 



ao72 = - (72)' - \ rfioi^O + i + i?3) - 1 i?2 + 1^ + 5^2 



^'72 



, (73 - 7l)/r for r ^ 



i 5,(73-71) 



9o7 



5oi?2 



for r = 

-727+^ ((73-7l)72 + (m-E3)) 

( - 2 (f] m + dioifi Oo) + / (1 - K (i^ </./2f ) 



+ (72 (-2 (i?l + i?3) - 3 m) - - (doR) ) 

7 ( - 2 (1 - K (Q 0/2)^) (m - ^)/r 
Kn{(f>i/r) {flo<P + <f>o^)) 



7 ( - 2 (^1 - K (f7 (/)/2)^) dr{m - M) 



+ K 



for r ^ 



for r = 



9ofi = fin 



agfio = 72f2o- 2 (73-71)% 

- 1 (^(m + B) + ^ ^2) (i-K (o<^/2)^))) 
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7(%/r + f Kf73</)(0i/r)) for r 7^0 



for r = 



1 



where m stands for 



9o% = -^[i~Kin (t)/2fj {Ri-m)n 

+ (^Mi- ^ (01 - ^) 0) - ^ (73 - 7I) f^o (B17) 

ao<^ = 00 (B18) 

ao0o = - 2720o-^(73-7l)-0f + 01 + 03 + % (B19) 

27 (0j_/r) for r 7^ 

V = <( (B20) 
— 2 7 9r0]_ for r — 

9001 = 01 - 03 - ^ (73 - 7I) 

Sorfioi- = - 3 72dioi-+ (m+ J0'f^of^dm-+5<iioia) (f^'^/2f) (B21) 
7f (0i/r)(17o0 + r!0o) for r 7^0 

7f (f^o + ^2 0o) 5^01 forr = 



= |(0r!i + 0ir!)(0i-03) 

I (0f}o + 00^^) (01 + 03) + J 01 (41^0 01 - 0f^ (S - ra) 
K 00 - 2 72 r^o + (^101-^ + (i?i + i?3 - ii?2) - ^ 

- i72r!0o-%0i J 



a2 



^ ^ ^ - % (m - :r3) + f^o (^(m + :R3-i^)--^i? 



2. The Constraints 



From (y) the foUowing set of constraint equations can be derived: 

5ieo° = -^(73-71) (B22) 

die = e — for r ^ 0, otherwise (B23) 

- r 

5i72 - i (ra - iTl) + J (^72 - ^ (7I + 73)) 

for r 7^ 0, otherwise (B24) 
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dvy = 7 



eo° + 7 



r 



r I R R2 \ 

+ 2 (-^ + i?i + ^ + i?3-dioiaf]- (71 + 73)721 

for r ^ 0, otherwise (B25) 
diR2 = 2j2 (^m-mj (diR) 

Sfl2 = 7 ( - 2 (1 - in<j>/2f) {Ri + m-^ m)/r 

for r ^ 0, otherwise (B26) 
6>iO = ill (B27) 



+ - 12'^ 1^2 00 + <^ (03 - (/.I) + Y (^^1 - ^3) J (B28) 

5i% = (^^-72)i7o-J%-§(|iri-:^+ir3j 

+ Jfl^ 1^201^ + ||03-01-720o-^0i + |(i?i-:^+ir3-g)jj 
for r 7^ (B29) 

= 01 (B30) 

5i0o = 01-03 + ^^^^01 (B31) 

■ — ^ ^ — ^ ^3 + 'yI R 
^101 = 01 + 03 + ^^-^00 + — (B32) 

aidioi- = ("ic + J 0' % 17 dioi- + SdHua) / (1 - « (O 0/2)') (B33) 
Sdioia = 7 (^3 7 ( (1 - « (0 0/2)') dioi^/r 

-^0i/r (fi0i+%0))^ 
for r 7^ (B34) 

where mc stands for 

mc = |(017i+0il7)(01 + 03)-|(017o + 0on)(0i-03) 



+ ^ 00 (4 f2g0i - 4 Oi 00 - 6 72 Oi + O (iJl - i?3)) 
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K 



l2 



On J 



(B35) 



Equation (Q), which must be fulfilled at one point, becomes 

Rfi^ + 6 mn^ -24: j2nno + 12^0^ -24- nrii- uQi^ - -rq'' 

- - r - - 4 

3 '7 

— K i?2 17^ 0^ - 6 K 72 ll"* 00 - 3 K fl"* 00^ - 6 K - (/) 01 
2 _ _ J, _ 

+3KrjVi^ = 0. (B36) 

Because of spherical symmetry there is also the identity 

nd,oi^^^^^ + {j2f-lR2-^. (B37) 



3. The regularity conditions 

At the center the following must hold for a regular spacetime: 

e = -7 = eii (B38a) 

73 = 7I (B38b) 

72=^(73 + 71) (B38c) 

Ri = m (B38d) 



i?2 = 2 i^Rl + R3j (B38e) 

dioi^ = (B38f ) 

= (B38g) 

01 = (B38h) 

01 = 03. (B38i) 

^0-^1=0 (B38j) 

M + /r3 = i- - d,o,n, - 1 8,n, + ( - 72' 

+ ^(iri-ir3)-o4+%(2:^72-(Jf 

+ ^.(5i^+(^-72) J)). (B38k) 
Equation (|B38j|) also guarantees that equation (B36) is fulfilled on at least one point in M, namely J^. 
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